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Abstract. The topological state of a two-dimensional triplet superconductor may be 
changed by an appropriate addition of magnetic impurities. A ferromagnetic magnetic 
chain at the surface of a superconductor with spin-orbit coupling may eliminate the 
edge states of a finite system giving rise to localized zero modes at the edges of the 
chain. The coexistence/competition between the two types of zero modes is considered. 
The reduction of the system to an effective Id system gives partial information on the 
topological properties but the study of the two sets of zero modes requires a two- 
dimensional treatment. Increasing the impurity density from a magnetic chain to 
magnetic islands leads to a finite Chern number. At half-filling small concentrations 
are enough to induce chiral modes. 
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Topological systems are the focus of current great interest. Topological superconductors 
oca have been studied in great detail, with particular emphasis on the existence of 
Majorana fermions at its edges [3j, or located around local perturbations, such as in the 
interior of vortices. A standard proposal is the one-dimensional Kitaev model of spinless 
fermions with p-wave pairing |lj, which displays topological phase in some parameter 
regimes, with the existence of localized zero energy modes, if the chain is f ini te (and 
long enough). Several two-dimensional superconductors also have topological properties 
such as the so-called p + ip pairing [SI El- Adding spin-orbit interaction and a magnetic 
held either through its Zeeman effect or due to the presence of vortices, a great variety 
of topological phases are predicted. 

Even though there are various candidates for a triplet superconductor, the proposal 
[2] of a conventional, very abundant in nature, s-wave singlet superconductor in 
proximity to a semiconductor wire in the presence of a Zeeman held [Hj, that also has 
topological properties, considerably increased the activity on this held. Experimental 
verihcations of the edge states have proved difficult but some experiments seem to 
provide good evidence for their existence, in particular a nanowire on top of a 
superconductor [9j and, more recently, a set of (classical) magnetic impurities on top of 
a conventional superconductor with their spin orientations arranged in some helical way 
HD], or a ferromagnetic chain in the additional presence of spin-orbit coupling m- Some 
controversy has however been raised [12., [THl O 15j, even though several predictions 
suggest their existence through different signatures, such as Andreev reflection [16| lTj. 

Since both a triplet superconductor and a magnetic chain induce topological states, 
we explore here the combined effect of the two by considering a set of magnetic impurities 
on top of a triplet superconductor. 

The model considered here has the same structure of that considered in [llj, with 
the difference that the underlying superconductor has spin triplet pairing. However, 
the inclusion of a Rashba like spin-orbit coupling implies the possible coexistence with 
spin-singlet pairing [18j, since parity is no longer conserved. If it exists it is assumed to 
be smaller and is neglected here. 

The system without the magnetic chain is described by the Hamiltonian [19j 



(1) 



where 


( 2 ) 


Hq — CfeCTo + H r . 


Here, e*, = — 2f (cos + cos k y ) — p is the kinetic part, t denotes the hopping parameter 
set in the following as the energy scale (t — 1), p is the chemical potential, k is a wave 
vector in the xy plane, and we have taken the lattice constant to be unity. The Rashba 
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spin-orbit term is written as 


Hr = s ■ cr = a (sin k y a x — sin k x a y ) , 


(3) 


where a is measured in the same units. The matrices cr x ,a y ,a z are the Pauli matrices 
acting on the spin sector, and <To is the 2x2 identity. 

The pairing matrix reads 




The spin S magnetic impurities act like local magnetic fields and are distributed 
along a chain at the center of the two-dimensional system, now considered finite with 
dimensions N x x N y . Their contribution to the Hamiltonian is written as 



(5) 
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where t- l are the Pauli matrices (i = x , y, z ) acting on particle-hole space. Here r 7 are 
the locations of the N < N x magnetic impurities distributed along the x direction as 

r . - ^Y»p 

J * 

2. Reduction to magnetic chain 

Since the main focus is on the low energy states of the system, we expect that these 
states will be zero energy modes, that will either appear localized at the edges of the 
magnetic chain or along the edges of the two-dimensional system, if open boundary 
onditions are used. As we will see, the choice of boundary conditions naturally affects 
the solutions found and it is the purpose of this work to determine the coexistence or 
competition between the low energy states. For the moment let us focus on the states 
associated more directly with the magnetic chain. This can be done by focusing on 
the degrees of freedom at the locations of the spin impurities. One possible way is to 
integrate the degrees of freedom of the rest of the system and to obtain an effective 
Hamiltonian for the remaining degrees of freedom. This can be achieved, for instance, 
in a path integral formalism and integrating out the degrees of freedom surrounding 
the magnetic chain, or using a lattice Green’s function method. A simpler method 
can be carried out by restricting the solutions of the BdG equations to the low-energy 
(subgap) states of the magnetic chain. As shown before for helical chains [20] and also 
very recently for ferromagnetic chains plus spin-orbit coupling [21 [ [22], the effective one¬ 
dimensional Hamiltonian will be of the type of the Kitaev model and naturally leads to 
zero energy edge states. 

We start from the Bogoliubov-de Gennes equations for the wave functions that in 
the momentum representation are written as 



( 6 ) 
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The 4-component spinor can be written as 

( u n {k ,t) \ 

Ur, \ Un(^kj .].) 

V n (~k, t) 
y v n (-k,i) 


(7) 


The impurity term is then an additional term to the Hamiltonian and will change 
the BdG equations. Introducing the Fourier components of the local operators at the 
impurity sites, we can write the BdG equations, including the magnetic impurities, as 


(w k - H k ) = - JS ^2 ( r z ® a z) Wj) e lkr F 

i 

The left hand side may be written as the inverse of the Green function 

G k 1 (w k )il> k = -JS ^2 ( T z ® cr z ) i/j(rj)e~ lkrj . 
j 


This leads to 


( 8 ) 

(9) 


'ijjk = — JS ^2 Gk(wk) [t z 0 a z ) ijj{rj)e lkrj . (10) 

j 

Taking the inverse Fourier transform we may get that the wave functions for the impurity 
locations are given by 

#G) = -JSj2 7^ elk{ri ~ ri)G k(w k ) (t z 0 a z ) ^(g)- ( n ) 

j * k 

Separating the j = i term from the right hand side we get that 

Fi/’(r i ) = -'^2,J(i,j)il>(r j ) (12) 

j#* 


where 


F = I + JS 


1 

Tv 


5 2 G k{w k ) ( t z 0 a z ) 

k 


(13) 


and 

J(i,j) = JS^=S2e ik ^-^G k (w k ) (t z 0 a,). (14) 

V N k 

Instead of attempting a direct solution of these equations we will follow previous 
treatments [20] that are suitable to obtain a solution for the low energy subgap states. 
Indeed, since the interest is in identifying the zero energy modes, a low energy expansion 
will be carried out. Also, we consider the limit of weak coupling to the impurities. This 
implies that in the term that connects the impurities, J(i,j ), we take the energy to be 
zero, e k = 0, and in the term F we take the limit of small energies, considering only 
a linear term in the energy. In this way we may get an effective low energy eigenvalue 
equation where the energy is the eigenvalue and the wave functions at the impurity 
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locations are the eigenfunctions. The point now is to identify the effective Hamiltonian 
for the magnetic chain that replaces the full 2 d original problem. 

In general, and recalling that the Green function is for the system without the 
magnetic impurities, we see that 

G k l {w) =w(tq® (Jo) - e k (t z ® cr 0 ) - as y (r 0 ® a x ) - ias x (t z 0 ia t ) 

+ d x (t x <g> a z ) - id y ( ir y ® (j 0 ) - A s (ir y ® ia y ) (15) 


Here s x = sin (Ay), s y = sin (Ay) and we are considering systems with no d z component. 
Let us now consider the spin singlet and the spin triplet cases separately. 

In the case of the spin singlet (d x = d y = 0, A s 7 ^ 0) we get that the inverse of the 
Green’s function is given by (in the limit of small energies, w) 


G = 


D 


e_ 

—ias + z(— 1 , — 1 ) 

2i Oi AgC^Sj^. 

—A s z(l, 1 ) 


ias-z(— 1 , — 1 ) 
e_ 

A s z(l, 1 ) 
i2aA s e k S- 


—2iaA s e k s_ 

A s z(l,l) 

e+ 

ias-z(— 1 , — 1 ) 


~A s z(l, 1 ) 
-i2aA s e k s + 
—ias + z(— 1 , — 1 ) 
e + 


\ 


(16) 


where e± = ±ekz(l,—l) — wz(l, 1) with z(r],rf) = e 2 k + gA 2 s -\-g'a 2 (s 2 x +s y ), s± = s x ±is y 
and D = z( 1, l) 2 — 4a 2 e^(s 2 + s 2 ). Note that the spin-orbit coupling (a ^ 0) induces 
elements of the spin triplet type (elements (1, 3), (2,4), (3,1), (4, 2)). 

Multiplication by (r 2 ® a z ) just changes the sign of all elements of the second and 
third columns. In the approximation of low energy subgap states, the interaction term 
between the impurities J(i,j ) does not depend on the energy, w, since we calculate 
the Green’s function at zero energy. We can rewrite the equation between the wave 
functions at the magnetic impurity locations as 


/ F J( 1,2) J(l, 3) ••• J(l, N) \ 


V>(1) \ 

J( 2,1) F J{ 2,3) ••• f(2,N) 


m 

J( 3,1) J( 3,2) F ••• J(3, N) 


V>( 3 ) 

\J(N, 1) J (N, 2) J(N, 3) ••• F J 


4>(N) ) 


= wJS 


1 

Vn 


E 


W-(l,1) 

D 


(‘ T z ® O z ) 


^(1) \ 

m 

■0(3) 


m) j 


(17) 


We see that the matrix on the left hand side acts like an effective Hamiltonian for 
the magnetic chain wave functions. As shown before, the presence of the spin-orbit 
interaction induces a pairing of the type d x , d y and the s-wave singlet pairing gets 
renormalized, as all the terms of the original Hamiltonian [23]. 

Considering now the case of triplet pairing and inserting the magnetic impurities, 
we have once again to find the Green’s function. This has now a simpler form, assuming 
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the absence of the spin-orbit interaction, since it is not needed to yield triplet pairing. 
It is enough to consider the structure of the Green’s function. This is given by 

/ — 6k — w 0 gL 0 \ 

0 —Cfe — w 0 — d + 

d .|_ 0 ek — w 0 

y 0 — gL 0 efc — w / 




(18) 


where here D = + d 2 x + d? y and d± = d x ± id y . We see that the structure of the triplet 

pairing is maintained and all terms get renormalized. 

From these results we see that the two-dimensional problem has been approximately 
reduced to a magnetic chain, with some effective Hamiltonian, that has the structure 
of a one-dimensional triplet superconductor, and therefore similar to the Kitaev model. 
Therefore one expects some regions in the phase diagram where the system may become 
topological and with the presence of edge states, if the magnetic chain is finite. However, 
as noted before, there are more complex terms with long range character both in the 
kinetic energy and in the pairing terms. These typically decrease with distance but may 
have a relevant effect. 


3. Phase diagram of magnetic chain 

The effect of longer range interactions has been studied before in the context of the 
Kitaev model [22 Eg m\ and in other topological superconducting systems mm- 
The Kitaev model is mapped by a Jordan-Wigner transformation to a spin problem, 
specifically the Ising model in a transverse field keeping only nearest-neighbor couplings, 
but it maps to more general spin interactions considering longer range terms. In these 
mappings the hopping and the pairing amplitude are related, but one may consider 
a more general case where they can vary independently. Taking, for instance, nearest 
neighbor and next-nearest neighbor terms (both in the hopping term and in the pairing), 
leads to a rich phase diagram that has been studied before. Interestingly, considering 
next-nearest neighbor terms leads to the possibility of two Majorana fermions at each 
end of the chain, and including longer range terms leads to an increasing number of 
Majoranas at the edges. In the limit where all first terms are absent and a given higher 
order term is the first non-vanishing, the system is equivalent to a set of copies of the 
original Kitaev chain (one for each decoupled sublattice) and therefore one has at each 
edge as many Majoranas as the number of copies of the system. The system of spinless 
fermions is in the BDI class with a Z number of Majoranas at the edges of the system 
if an effective time reversal symmetry holds and the states are protected due to the 
existence of a chiral symmetry, which in this problem is just the particle-hole symmetry 
of a BdG problem. 

This has been determined before for the Kitaev model composed of spinless 
fermions. The case of the magnetic chain is similar, but not quite the same. In this 
problem we have spinful electrons and interactions that are more complex than just 
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{t=l,d=0.6} {t=l,d=0} {t=0,d=0.6j {t=0,d=0} 
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Figure 1 . (Color online) Phase diagrams for the extended Kitaev chain with second 
neighbor hopping and second neighbor pairing, indexed by the winding number. 


the triplet pairing of the Kitaev model. We will determine the phase diagram of these 
magnetic chains and also present some results, for the sake of completeness, for the 
Kitaev chain. A chiral symmetry is also present enabling to classify the topological 
phases by a winding number. 

A direct method to determine the regions in the phase diagram that are topological 
and how many Majoranas are obtained at each edge, is to diagonalize the Hamiltonian 
and determine its eigenstates and eigenvalues looking for zero modes. Another, more 
systematic way, is obtained calculating the winding numbers that give directly the 
number of edge states in the problem [23 m\. This is somewhat more reliable since 
the diagonalization is perfomed for a finite system and, even though the energies of the 
Majorana states are rather low (in some cases of the order of 10 -16 ), in other cases they 
are higher (of the order of 10 -5 ), but still the winding number is nonvanishing. 

3.1. Id extended Kitaev model 

We start with the case of the extended Kitaev model. The model in real space may be 
written as 

H = -1 (aW+ 4+i a) (4a+ 2 + 4+2 a) - i-i 4 a 

iii 
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+ a ^ (ciC i+ 1 + 4+14) + ^2 ( c * c *+ 2 + 4+24) (is) 

i i 

where t and t! are the nearest and next-nearest neighbor hoppings, respectively, and A 
and A' the nearest and next-nearest neighbor pairings, respectively. This Hamiltonian 
may be diagonalized for a finite system of size N x using open boundary conditions, and 
yields the presence of edge states of zero energy. 

In momentum space the model is simply written as 


where 


& = \ £ ( 4 . c-k) 

k 



H k 


£fc - h 

—i A sin A: — i A' sin 2k 


i A sin k + i A' sin 2k 
~£k + P 


( 20 ) 


( 21 ) 


with Cfc = —2tcosk — 2t'cos2k. 

Straightforward diagonalization gives = (—+ p) 2 + (A sin k + A'sin2A:) 2 . 
Transitions between different topological phases occur when the system becomes gapless 
which implies — e*, + /i = 0 and A sin k + A' sin 2k = 0. This occurs when cos k = 
—A/(2A'), which leads to t' = — A(pA' — tA)/(A 2 — 2(A') 2 ). If t! = 0, A' = 0, there are 
zeros for k = 0 and p = —21 or k = n and p = 21. If only A' = 0 we get zeros at k = 0 
and 2t + 2t! + p = 0 or k = 7r and —2 1 + 2t' + p = 0. If only t' = 0, we get p = tA/A'. 

To calculate the winding number we follow the usual way by identifying a matrix 
that anti-commutes with the Hamiltonian matrix. This is obtained verifying that 
r x HkT x = —H. Defining a matrix, T, composed by the eigenvectors of t x , we can 
reduce the Hamiltonian to an off-diagonal form 

™ r ’ = (m q V ) (22 > 

where, in this case, q{k) is just a function given by 

q{k) = efc — p — i A sin k — iA' sin 2k. (23) 


The winding number may be defined by 

1 = h L Tt - W )' 1 J+ f m) • (24) 

The phase diagram where the various topological phases are characterized by the 
winding number is represented in Fig. [lj where the results are presented for various 
values of the chemical potential p = 0,1, 3. In the case of only nearest-neighbor hopping 
and pairing, and p = 0,1, we are in topological phases and if p. = 3 we are in the 
topologically trivial phase (we take t = 1, d = 0.6). Also, we consider cases in which one 
or both of the nearest-neighbor terms are absent, namely t = 1 ,d = 0, t = 0 ,d = 0.6 
and t = 0, d = 0. 


Magnetic chains on a triplet superconductor 


9 



h z /2,- 

Figure 2. (Color online) Phase diagram for an effective magnetic chain indexed by the 
winding number for a s-wave superconductor (left panel) and a p-wave superconductor 
(right panel). 


The effective interaction J(i,j ) decays with distance and, as a result, the nearest- 
neighbor terms are expected to be smaller than the nearest-neighbor ones (t'/t < 
1,A'/A < 1). The results are however presented beyond this region for completeness. 
In the limits of high values of t'/t and/or A'/A one expects changes with respect to the 
original Kitaev model, as discussed above. However, the results in Fig. [l] show that, 
even in the regime where t'/t< 1, A'/A < 1, the winding number changes with respect 
to the Kitaev model, and the topology and number of edge states changes. For large 
values of \t'\ or \d'\, the winding number is W = 2, —2 indicating two edge states (on 
each end of the chain). Around the central region W = 1, 0, —1, as expected of the limit 
t'/t « 1, A'/A << 1. As shown here and shown before by other authors, the topology 
is changed in some cases even if the additional terms are small, rendering the system 
trivial or topological. If t = 0 the system is either trivial (W = 0) or has two edge states 
(W = 2, —2) for the various chemical potential values. As the chemical potential grows 
the trivial phase becomes more widespread. 


3.2. Id magnetic chain 


We consider now a magnetic chain described by an Hamiltonian matrix of the form 

( — h z iasink —idsink A s \ 

—iasink Ck + h z — A s —id sink 
id sink —A s —Ck + h z —iasink 


H, = 


(25) 


y A s id sink iasink —ek~h z J 

where e*, = — 2t cos k — /i and, as above for the two-dimensional superconductor with 
magnetic impurities, we consider the presence of spin-orbit coupling a, and triplet and 
singlet pairings along the direction of the chain; h z is a local magnetic field that simulates 
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the ferromagnetic ordering of the impurities. This model is just the one-dimensional 


reduction of the 2d Hamiltonian considered in Eq. [T] (taking k = k x ). 


The operator that anti-commutes with the Hamiltonian is now the matrix t x £3 (To, 
and we obtain that q{k ) is a 2 x 2 matrix given by 



e k ~ h z + id sin k ia sin k — A s 

—ia sin k + A s £& + h z + id sin k 


(26) 


where here = —2 1 cos k — ji. The winding number may be calculated in the same way. 

The results are shown in Fig. | as a function of y and h z , for a s- 
wave superconductor and for a p-wave superconductor. In the case of the s-wave 
superconductor the winding number vanishes if the spin-orbit coupling is absent. In the 
p-wave superconductor, even if the spin-orbit coupling vanishes, the winding number is 
non-vanishing in regions of the phase diagram. If a ^ 0 then the phase diagram does 
not depend on a, for both pairings. 

The phase diagram of the triplet chain is somewhat similar to the triplet 2d system 
[19J with a rescaling of both the chemical potential and the magnetic field from 4 1 to 
2t. In the 2d case and p > 0, there are two topologically distinct phases with (7 = 0 
and (7 = —2, where (7 is the Chern number, but both with two edge states. The same 
occurs in the magnetic chain. For y > 2 (and small h z ) and h z > 2 (small /j) , there are 
no edge states and the system is trivial. There is also a region with W = — 1 (one edge 
state) as for the 2d case. 

The s-wave superconductor is however different from its 2d counterpart, since the 
phase with two edge states (IT = ±2) does not appear and only the phase with IT = — 1 
is topologically non-trivial. 

4. Two-dimensional system 

We consider now a finite two-dimensional system of dimensions N x x N y , along a 
longitudinal direction :c, and a transversal direction y. We apply different types of 
boundary conditions along the longitudinal and transverse directions leading to different 
sets of edge states. We write 



(27) 


and rewrite the Hamiltonian matrix in terms of the operators (27) as 



Jx Jy 


''Pjxdy ,t \ 


r jx -jy 4 

jxjy „/,t 


( 28 ) 
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The operator H lr Jy reads 


where 


H- ■ = \ A B 

12 Jxjy \ Q jj 


A = 


B = 


( —h z — ep — tg+ — trf + | r t ft + 


C = 


D = 


-f rj^ + YiV- h z — ep — trf^. — trf + 

0 


ryyC d /yyJ 

' 2 'I- 2 i 'I- 

0 


_d x I d-^y 

2 'I- b 2 i 'I- 


d x _ d_ rn y 

— I 2'I- 2 i'l- 

q ct ^x d ^y 



h z + e F + tri x + + trf + 

— | xy 2 

2 'I- ' 2z ' — 


a^x \ a^^y _^ 


_ot L X I 

2 , /- ' 2z 'I- 
+ + £??+ 


where • U« ~ '-'.l • ij/ and ^L jy V±^j x j y = 'ij/Tj.-i + 

The diagonalization of this Hamiltonian involves the solution of a 
(47V,Ai y ) x (47V,7V y ) eigenvalue problem. The energy states include, in general, states 
in the bulk and states along the edges. The local magnetic field describes different 
distributions of magnetic impurities, as needed. 


— ft! ,t 


d 


t 


(29) 

(30) 

(31) 

(32) 

(33) 


— „/,t 


4-1. One impurity 

We begin by considering the effect of a single magnetic impurity on the p- wave 
superconductor and compare the results with a magnetic impurity on a s-wave 
superconductor. 

The effect of a single impurity in a conventional s-wave superconductor is well 
understood EH- Since the impurity spin acts like a local magnetic field, the electronic 
spin density will align along the local spin. For small values of the coupling there is a 
negative spin density around the impurity site. At the impurity site it is positive, as 
expected. For larger couplings the spin density in the vicinity of the impurity site is 
positive. At small couplings the many-body system screens the effect induced by the 
impurity, inducing fluctuations that compensate the effect of the local field, in a way 
that the overall magnetization vanishes. However, for a strong enough coupling the 
many-body system becomes magnetized in a discontinuous fashion. One interpretation 
is that, if J is strong enough, the impurity breaks a Cooper pair and captures one 
of the electrons, leaving the other electron unpaired, and thus the overall electronic 
system becomes polarized. The impurity induces a pair of bound states inside the 
superconducting energy gap, one at positive energy (with respect to the chemical 
potential), and another at a symmetric negative energy. Considering a higher value 
for the coupling one finds that the levels inside the gap approach the Fermi level. There 
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Figure 3. (Color online) Lowest energy value as a function of the coupling between 
the impurity spin and the spin density of the electrons. 


is a critical value of the coupling for which the two levels cross in a discontinuous way 
such that it coincides with the emergence of a finite overall magnetization. After the 
level crossing occurs, the nature of the states changes and, as the level crossing takes 
place, the spin content also changes. The level crossing occurs between one state that 
describes an uncompensated local spin (at smaller coupling) and a state where the 
impurity spin is compensated (partially since for the classical description to be valid the 
spin has to be large). 

In Fig. [3] we show the results for the lowest (positive) energy state for the swave 
and the p -wave cases and for different chemical potentials // = 0. —3. In both cases, as 
a function of the coupling J, the lowest energy state approaches zero and then increases 
again, signalling the quantum phase transition (QPT). Note that here the calculation is 
not done self-consistently and the superconductivity is assumed constant, as the result 
of a proximity effect (A s , d = 0.6, respectively). In a full self-consistent solution the gap 
function is decreased at the impurity site and at the quantum phase transition changes 
sign ( 7 r shift) and the energy level does not vanish in the s-wave case. In the p-wave case 
the quantum critical point occurs at higher coupling values and for p = 0 it is rather 
flat up to high coupling values. Note that the effect of a single impurity is significant 
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Figure 4. (Color online) Lowest state wavefunctions for the case of OBC along x 
and y directions for g = 0, left column, and g = —3, right column, for h z = 2,3,5 
for a system of size 210 x 19 and a chain of magnetic impurities of length 150 sites 
(■ N<N X ). 


in the properties of the system as shown before [HU EU ESI HI] , including its effect on 
transport properties inducing, for instance, currents [35] and an anomalous Hall effect 
[36], if spin-orbit is present. 

4-2. Magnetic chains 

Increasing the number of impurities the number of in-gap states increases (two per 
impurity) and the gap gets filled. If the impurities are arranged in a chain and if this 
is long enough, Majorana edge states have been proposed to occur, if for instance the 
alignment is parallel and there is spin-orbit coupling HU- 
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Figure 5. (Color online) Lowest state wavefunctions for the case of OBC along x 
and PBC along y, for /j = 0. left column, and /j = —3, right column, for h z = 1, 2, 5 
for the same system of Fig. [4j 


It is therefore interesting to see what happens if the magnetic impurities are 
deposited on top of a p -wave superconductor. One expects two types of edge states. 
Due to the magnetic chain one expects localized edge states at the ends of the chain 
that, if the chain is long enough, are zero energy Majorana modes. On the other hand, 
a p -wave superconductor at zero magnetic held has propagating edge modes along the 
borders of a finite system in a stripe geometry (with periodic boundary conditions along 
one direction and open boundary conditions along the other). In order to observe the 
edge states along the chain we choose open boundary conditions along the x direction 
(magnetic chain direction) and along the y direction we choose either open boundary 
conditions or periodic boundary conditions. Also, we consider a case where the magnetic 
chain is imbedded inside the two-dimensional system (N < N x ) and another case where 
the magnetic chain extends up to the edges of the system (N = N x ). In the first case 
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a spatial separation between the two types of low energy edge states is possible while 
in the second case, if both exist, they will be superimposed. In both cases there is a 
competition between the two types of states, being more evident for the second case. 

We consider first the case where the magnetic chain extends from site 30 up to site 
150 along the x direction, centered in the middle of the system in the y direction, in a 
system of dimensions 210 x 9. The lowest energy wave functions are either symmetric 
in x or, if not, there is another degenerate state or a corresponding state at negative 
energy that is its mirror. Since both types of edge states are rather localized, for better 
visualization we focus on a region close to one border of the two-dimensional system. 
Some states have very small energies but other edge states have subgap energies large 
but still clearly smaller than the states in the continuum, typically at least one order of 
magnitude smaller. 

It is perhaps helpful to recall the phase diagram for a two-dimensional system with a 
uniform magnetic field mm- In this work we have a non-homogeneous distribution of 
magnetic fields located at the magnetic impurities, since these interact with the system 
as a local Zeeman term. At pi — 0, —3 and zero magnetic field, the system is in a Z 2 
phase with two degenerate edge modes on each edge, considering a stripe geometry. 
Turning on a magnetic field the system at y, = 0 changes to a Z phase with Chern 
number C = —2, for any magnetic field between 0 < h z < 4. At h z = 4 there is a 
quantum phase transition to a trivial regime with (7 = 0 and no edge states. In the 
case /i = — 3 if 0 < /^ < 1 the system has two edge states but (7 = 0. In the region 
1 < h z < 3 the system has one edge state and (7 = 1, between 3 < h z < 7 again 
one edge state but with (7 = — 1, and above that magnetic field the system becomes 
trivial with (7 = 0 and no edge states. In every change of the Chern number there is a 
quantum phase transition and the bulk gap closes. The number of edge states may also 
be understood calculating the topological entanglement entropy |38j 

In Fig. [4] we apply open boundary conditions (OBC) along both spatial directions 
and show the density of the lowest state eigenfunction, as a function of x and y, 
for two values of the chemical potential y = 0, —3, and consider that each magnetic 
impurity can be seen as a local magnetic field with values h z = 2, 3, 5. Also, we take 
a = 0, d = 0.6, A s = 0. Consider first y = 0. In the cases of h z = 2,3,5 the groundstate 
is four times degenerate, with an energy of the order of 10 -7 or 10 —6 , and the next 
set of states has an energy of the order of 10 -2 . In the left column of Fig. [ 4 ] one of 
the degenerate states is shown. In this case both types of edge states are clearly seen. 
The state along the border of the system peaks at the edges along y and the state 
associated with the end of the magnetic chain is peaked at this edge and decays towards 
the edges along the y direction. Both states are very localized along the x direction. If 
the chemical potential is y = —3, the border state dominates and any peak at the edge 
of the magnetic chain is hardly visible at low magnetic fields, but its weight increases 
with the magnetic field being clearly visible at h z = 5. However, in this case the two 
peaks overlap and in the region between them the wave function is clearly finite and 
appreciable. The lowest energy state is non-degenerate and with very small energy of 
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Figure 6. (Color online) Lowest state wavefunctions for the case of OBC,OBC for 
p = 0, left column, and p = —3, right column, for h z = 2,4,5 for the edges of the 
magnetic chain at the border of the system (TV = N x ). 


the order of 10“ 15 Actually, one may note that the two peaks somewhat always overlap 
through the borders of the system along the y direction [39]. Increasing the size along 
y does not diminish the spatial extent along y. 

Changing the boundary conditions along y to periodic (PBC), as shown in Fig. 
[5j tends to separate the two types of edge states. In the case of // = 0 the peak 
asociated with the magnetic chain dominates (with now some overlap with the border 
state extending along x ) and in the case of p = —3 the border states are even sharper. 
The border states are clearly invariant along y due to the PBC but the edge states at 
the chain ends are still peaked at the middle of the system. Note, however, that in 
the case of p = 0 the lowest energies are now quite high. The groundstate is doubly 
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degenerate with an energy of the order of 10 -2 . In the case of /.i = —3 the groundstate 
is now doubly degenerate with a rather small energy of the order of 10 -14 for h z = 1, 2 
and non-degenerate with an energy of the order of 10 -18 at h z = 5. 

Extending the magnetic chain up to the border of the system reverses the results 
between the two chemical potentials: in the case of ji = 0 the border state dominates, 
while if n = —3 the chain edge state dominates, as evidenced by the peak along y 
centered at the chain location. These results are shown in Fig. [6} Again for // = 0 the 
energies of the lowest states are rather high, of the order of 10 -2 for h z = 2,4, and of 
the order of 10 -3 for h z = 5, while for /i = —3 they are nondegenerate for h z = 2,4 of 
the order of 10 _n or smaller, while for h z = 5 the energy is of the order of 10 -2 and the 
state is doubly degenerate. 

Associated with edge states we expect a bulk topological invariant. For 
translationally invariant systems a classification of the various possible symmetry classes 
has been extensively studied [TO m\. In the presence of magnetic field time reversal 
symmetry is broken and so a non-vanishing Chern number may be used to classify the 
phases. This has been used in the case of a uniform magnetic field, as reviewed above. 
In the case of the magnetic chain immersed in the two-dimensional space, one may 
wonder if the magnetic field applied on the system is enough to change its topology. 
The usual way to calculate the Chern number, obtained through an integration over the 
Brillouin zone, does not apply here since there is no translational invariance. However, 
as discussed in the following subsection, one may define a Chern number in real space 
by integrating over twisted boundary conditions. This method may be used for non- 
homogeneous situations where crystal momentum is not a good quantum number. A 
classification for systems with defects has also been carried out [42 j. 

4-3. Topological invariant: Chern number in real space 

A possible choice of a topological invariant to characterize each topological phase 
is obtained calculating the Chern number associated with integration over twisted 
boundary conditions in the two-dimensional system [43J. 

If a quantum system has an Hamiltonian that depends on some parameter R and 
if it is periodic in this vector, than it can be shown that the integral of the z component 
of the Berry curvature over the surface the vector sweeps, S, is an integer called the 


Chern number 

Cn = ^~ f dSn^(R). (34) 

2tt Js 

The index n respects to a given eigenstate of the adiabatic energy levels 

H(R)\u n (R)) = E n (R)\u n (R)). (35) 

The z component of the Berry curbature of such a level is defined as (here (1,2) = (x, y) ) 

n ' n| = wk A * )(R) ~ w 2 A 3 {R) <36) 




Magnetic chains on a triplet superconductor 


18 


where the Berry connection is defined as 

A^(R) = (u n (R)\iV R \u n (R)). (37) 

In a translationally invariant system the parameter R may be chosen as the momentum 
and the integration runs over the Brillouin zone. In the problem at hand, which is not 
translationally invariant, the procedure may be chosen in terms of twisted boundary 
conditions in both spatial directions. 

The integration over the Brillouin zone of the Berry curvature has in general 
numerical problems as one approaches a transition point (for which the gap between 
bands closes) and it is more convenient to calculate the Chern number in the lattice 


The twisted boundary conditions are taken as 
«40 + N,a,) = e a ‘i4,(r) 


(38) 


where 0 = ( 61 , 62 ) and H( 6 )\u e n ) = E^\u e n ), and cp are the basis vectors of the direct 
lattice, along direction i. The problem is solved for a finite system with size Ni x N 2 , 
typically taken large enough. In general there are degeneracies in the energy spectrum 
which render the calculation of the Chern number difficult and so a method on the 
lattice is also prefered. One considers the Slater determinant of the single-particle wave 
functions of the Hamiltonian H( 6 ). The groundstate may then be represented by a 
matrix 


II 

(tif 

rtf ■ 

■ 

■■ K e \ 
•• 



^ ■ 

•• ) 


(39) 


where N = NiN 2 and M is the number of points M = Mi AT in the space of the twisted 
boundary conditions. The lattice Chern number may then be obtained in a way similar 


Mi M2 

C = — ^2 arg ((^ ei \^ dl+1 )(^ dl+1 \^ dl+1+2 )(^f dl+1+2 \^ dl+2 )(^f dl+2 \^ dl )) (40) 
1=1 

where the vectors 1 and 2 denote the two spatial directions corresponding to the two 
types of twisted boundary conditions along the system’s spatial directions. Since each 
state is now a many-body state given by a Slater determinant, the overlaps between two 
states with different boundary conditions are given by 

(T^') = det (41) 

and the Chern number is obtained by 

1 M 

c = ^z ^2 argX p ( 42 ) 

71 1=1 

where X p are the eigenvalues of the matrix product 

M 

1=1 


( 43 ) 
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Figure 7. (Color online) Chern number as a function of concentration and magnetic 
field for magnetic islands of increasing size (concentration) for g = 0 and a = 0, 0.6 
(left and right panels, respectively). 



Figure 8. (Color online) Chern number as a function of concentration and magnetic 
field for magnetic islands of increasing size (concentration) for g = — 3 and a = 0, 0.6 
(left and right panels, respectively). 


However, this procedure is very time consuming. It has been shown that it is enough to 
calculate the Chern using only periodic boundary conditions and not twisted boundary 
conditions [45]. The Chern number is obtained by a similar expression but the 
eigenvalues X p are the eigenvalues of the M x M matrix 

F = C qoqi C qiq2 C q2q3 C q3qo (44) 


where the momenta are defined by 

z N , 2tt x , 2tt 2tt x , 2n x 

«° = I 0 ' 0 )’*'- = j ^);«3 = 

The matrices are defined by its matrix elements, mn, as 


(45) 


rrmn 

^q,q' 


E 


9 = 0 \ 


O i(q-q')-rij,n,e=o 


(46) 


Applying this method to the problem of the magnetic chain yields zero Chern 
number, unless the magnetic field is very large, since the magnetic impurity 
concentration is small, and if we can think in a mean-field like argument the average 
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Figure 9. (Color online) DOS for various concentrations c, = 

0.02,0.11,0.27,0.51,0.82,1 for g = 0 and various magnetic fields h z = 0,1,2,3,4,5. 


magnetic field is small. One expects therefore that in the case of pi = —3 the Chern 
number should vanish. In the case of p, = 0, since in a uniform magnetic field any small 
field will change the Chern number from (7 = 0 to (7 = — 2 the Chern number may be 
finite. In order to better understand the effect of the impurity concentration, we consider 
in the next subsection magnetic islands of growing size, centered in the 2d system. A 
system of various chains has been studied recently for a singlet superconductor [46]. 
Also, edge states at the borders of a superconducting island in a topological insulator 
[42] and a magnetic island on a singlet superconductor [48j have been studied very 
recently. 

4-4- Magnetic islands 

We consider squares of increasing size centered in the two-dimensional system with 
dimensions (2 m + 1) x {2m + 1), where m = 1, • • •, 10 in a system of 21 x 21 sites. 
The Chern number can then be calculated for these regular distributions of magnetic 
impurities. The results are shown in Figs. 01 as a function of the magnetic field due 
to each impurity and the concentration (q = (2 m + 1) 2 /21 2 ). 

The case of concentration q = I is the case of a uniform magnetic field. In this case 
the results using the real space representation of the Chern number, agree remarkably 
well with the results obtained in momentum space. Lowering the concentration the same 
trend is still observed. At zero magnetic field clearly (7 = 0. At small concentrations the 
calculation of the Chern is unstable and gives results that are between the results for 
q = 1. The case of fi = 0 is more stable, as expected. The unstable results for the Chern 
occur particularly in the regions where the Chern is changing its value. However, the 
results for the Chern are often half-integer indicating non-analyticies, as expected since 
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Figure 10. (Color online) DOS for various magnetic fields h z = 1,2,3,4,5 and 
various concentrations c t = 0.02,0.11,0.27,0.38,0.51,0.66,0.82,1 for g = 0. 


there is a quantum phase transition. In the uniform case this is a sharp transition with 
the closing of the bulk gap but, for the non-homogeneous case, the transition regions are 
quite extended, particularly as the concentration decreases, the magnetic field is large 
and particularly for // = —3 (as compared to p = 0). These instabilities are associated 
with states in the gap approaching zero, as will be seen in the calculations of the density 
of states. 

In Figs. |9pO|II1 we show results for the density of states (DOS) for various cases. 
We see that often when a gap closes (or nearly closes) the energy levels approach 
these low energies in pairs centered around zero energy. This splitting of the bands 
presumably gives rise to Chern values that are half-integer. Something similar occurs 
in the calculation of the Hall conductance of a vortex lattice at high magnetic fields 
in the Landau level description calculated in the diagonal approximation, where as the 
gap amplitude approaches zero the Emitting value of the Hall conductance does not 
approach the normal state (integer) result, but a half-integer value. This result is due 
to the doubling of the energy bands and associated filling of half the bands around the 
Fermi level [49J. Something similar happens in graphene [50j. 

In Fig. [9] we consider /i = 0, various concentrations and magnetic fields. The case 
Ci = 1 corresponds naturally to a uniform magnetic field. At h z = 4 the gap closes 
due to the QPT between C = — 2 and (7 = 0 and there is a zero energy peak in the 
DOS. Introducing a small concentration of magnetic impurities, states are introduced 
in the gap that move to lower energies as the magnetic field at each impurity location 
increases. Increasing the impurity concentration (the size of the magnetic island) the 
behavior is no longer monotonic. As the field increases, the subgap states move away 
from small energies and have values that exceed the gap at zero magnetic field. In the 





































































Magnetic chains on a triplet superconductor 


22 






Figure 11. (Color online) DOS for various magnetic fields h z = 0.5,1,2,3,4,5 and 
various concentrations c t = 0.02,0.11,0.27,0.38,0.51,0.66,0.82,1 for g = —3. 


uniform case (c t = 1) the energy gap is larger than the zero-field value. The presence of 
low-energy states at small concentrations leads to non-analyticities and a Chern number 
that is not well-defined, as shown in Fig. [7| At larger concentrations the subgap density 
of states moves from small energy to higher energies. The vicinity of h z = 4 is however 
close to the QPT at c t = 1, and signals the transition between C = — 2 and C = 0. even 
at small and intermediate concentrations. For h z > 4 the system is in a trivial phase 
independent of the impurity concentration, as clearly illustrated in Fig. [7j 

In Fig. [lO] we also consider /i = 0 but each panel has a fixed magnetic field and 
the impurity concentration is varied. The case h z = 4 shows a filled gap for all q, as 
discussed above; the same occurs for h z = 5, but, at least for high concentrations, the 
system is clearly in the trivial phase, as shown by the value of C = 0 of Fig. [7| Except 
at the transition points of the QPT, the gap in the uniform case is large. As h z increases 
the low-energy part of the gap gets increasingly filled but in a non-monotonic way as 
a function of the impurity concentration. At the small value of h z — 1 we see that, 
excluding the lowest concentrations, increasing the size of the magnetic island the low 
energy states increase in energy exceeding the value at zero field. The same trend occurs 
for higher fields, but low energy states are found as one enters the trivial phase h z > 4. 

Finally we consider // = —3 in Fig. [llj The results for the Chern number in Fig. 
[8] are qualitatively the same, but now there are two QPT for C — 0 to C = 1, and to 
C = — 1. At Cj = 1 the first transition occurs at h z = 1 and the second transition at 
h z = 3. The first transition is clearly seen for c, = 1 with a zero-energy peak. However, 
the transition at h z = 3 is not seen. In the thermodynamic limit the first transition 
occurs at k x = 0 while the second one at k x = tt m l5l] . The results shown are 
obtained for a system of size 21 x 21. The first transition is captured by such a small 
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system but the second one has the wrong parity. We have checked that a system of 
size 22 x 22 captures the second transition (as it also does capture the first one) with 
a very small energy of 10 -15 . Note the splitting of levels at low energies that explains 
the half-integer values for the Chern number. At small fields (h z < 1) the dependence 
on the impurity concentration is now the opposite of the case at half-filling: a small 
concentration induces states that are of an energy that is larger than the zero-field 
case and, as the concentration increases, the energies get smaller, but do not differ 
appreciably from the zero-field case. Crossing the first QPT states appear in the gap 
down to very small energies, again in a non-monotonic way. At and after the QPT at 
h z = 3 (in the uniform case) there is a high density of states near zero energy that 
remains in the trivial phase. This happens as long as the magnetic island does not fill 
the entire system. Close to q = 1 the gap increases fast and, in the uniform limit, it is 
large. Considering larger systems changes the details of the energy levels distributions 
but the qualitative analysis remains the same. 

5. Summary 

The change of topology of a triplet superconductor due to the addition of magnetic 
impurities was studied in this work. Under appropriate conditions the zero modes 
along the borders of the superconductor are replaced by edge states localized inside the 
superconductor or at the border and at the ends of a magnetic chain if long enough. 
Increasing the density of impurities to magnetic islands of growing dimension leads to 
chiral states that replace the helical states at zero magnetic field, as shown by the 
non-vanishing Chern number. 

The author acknowledges discussions with Eduardo Castro on the Chern number 
in real space and support from FCT through grant UID/CTM/04540/2013. Partial 
support from the Japan Society for the Promotion of Science is gratefully acknowledged, 
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